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ABSTRACT 

Context. According to time-distance hehoseismology, information about internal fluid motions is encoded in the travel times of solar 
waves. The inverse problem consists of inferring three-dimensional vector flows from a set of travel-time measurements. While only 
few tests of the inversions have been done, it is known that the retrieval of the small- amplitude vertical flow velocities is problematic. 
A thorough study of biases and noise has not been carried out in realistic conditions. 

Aims. Here we investigate the potential of time-distance hehoseismology to infer three-dimensional convective velocities in the 
near-surface layers of the Sun. We developed a new Subtractive Optimally Localised Averaging (SOLA) code suitable for pipeline 
pseudo- automatic processing. Compared to its predecessor, the code was improved by accounting for additional constraints in order 
to get the right answer within a given noise level. The main aim of this study is to validate results obtained by our inversion code. 
Methods. We simulate travel-time maps using a snapshot from a numerical simulation of solar convective flows, realistic Born travel- 
time sensitivity kernels, and a realistic model of travel-time noise. These synthetic travel times are inverted for flows and the results 
compared with the known input flow field. Additional constraints are implemented in the inversion: cross-talk minimization between 
flow components and spatial localization of inversion coefficients. 

Results. Using modes /, pi through 774, we show that horizontal convective flow velocities can be inferred without bias, at a signal-to- 
noise ratio greater than one in the top 3.5 Mm, provided that observations span at least four days. The vertical component of velocity 
(v^), if it were to be weak, is more difficult to infer and is seriously affected by cross-talk from horizontal velocity components. We 
emphasise that this cross-talk must be explicitly minimised in order to retrieve in the top 1 Mm. We also show that statistical 
averaging over many different areas of the Sun allows for reliably measuring of average properties of all three flow components in the 
top 5.5 Mm of the convection zone. 
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1. Introduction 

The sub-surface of the Sun is optically thick, preventing us 
from directly observing the interior layers. Understanding the 
properties of the plasma in these regions has consequences for 
the theories of convection, stability of sunspots, the dynamics 
of stratified convection, and others. Most current knowledge 
about conv ection comes primarily from computati onal work 
(e.g. Vogler et al.| [2005[ [Benson et all [^06; Re mpel et aL 
|2009| ). Helioseismic inversions of the sub- surface flows will play 
an important role in constraining these theories. 

A powerful way of imaging the solar interior is via infer- 
ences gathered from studying the statistics of the acoustic and 
surface gravity waves at the surface. Solar pressure and surface 
gravity modes are generated randomly by the vigorous turbu- 
lence in the upper convection zone. These oscillations are ob- 
served in the solar photosphere by measuring Doppler shifts 
of photospheric absorption lines. Forward modelling allows us 
to relate anomalies (like flows, thermal hot/cold spots etc.) to 
changes in helioseismic observables. 
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The aim of helioseismic inversions is to reveal the structure 
of the subsurface flows (rotation, meridional circulation, convec- 
tion), magnetic fields, and to measure deviations in the plasma 
state parameters (temperature, density, pressure) from a quiet 
Sun average. In this paper, we focus on travel times ( [Duvall 
ft al.nl993| ), i.e., quantities that emerge from fits to cross cor- 
relations of observed signals. Time-distance hehoseismology is 
used to measure and interpret changes in travel times of seis- 
mic waves caused by inhomogeneities in the structure of the 
Sun (see review by [Gizon et al.[|2QTQ| ). In recent years, time- 
distance he hoseismolog y has been us ed to invert for near- surface 
flows (e.g., [Gizon et al., 2000 ; Duvall & Gizon[ [20001 piao^ 
|Kosovichev| [2004t [Jackiewicz et al.| [2008|), for flows beneath 
sunspots (e.g., [Duvall et al.[ [1996| [Zhao et al.|[20 0H 'Couvidat 
eFaL] [20061 [Cameron et aLll2Q08llGizon et al.[[2 009; Moradi 



et al. [201 0) and flows in their vicinity ^izon et aL ^2000), study 



the rotational gradient at t he base of the convection zone (e.g., 
Hanasoge & Duvallj [20091 ), etc. 



Helioseismic inversions are performed using two principal 
methods: The regularised least squares (RLS) and optimally lo- 
calised averaging (OLA). The RLS met hod (in time-dista nce he- 
hoseismology used for the first time by Kosovichev 1996j ) seeks 
to find the models of the solar interior, which provide the best 
least-squares fit to the measured travel-time maps, while reg- 
ularising the solution (e.g., by requiring the smooth solution). 
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The OLA method wa s developed for geoseismology d Backus & 
|Gilbert| |1968[ |1970| ). A form suitable for use in helioseismol- 
ogy was devised by |Pijpers & Thompson ( 199 2), who formu- 
lated the Subtractive-OLA method. SOLA is based on explicitly 
constructed spatially confined averaging kernels by taking linear 
combination of sensitivity kernels, while simultaneously keep- 
ing the error magnification small. The resulting coefficients are 
then used to linearly combine the travel-time maps and obtain 
the estimate for structure and magnitude of solar plasma per- 
turbations. A SOLA-type inversion is the principal method dis- 
cussed in the current paper. The SOLA has been used in time- 
distance local helioseismology in the past by J ackiewicz et al. 
pQOTf |20Q8 ) who demonstrated the ability of SOLA inversions 
to reveal the structure of the 3-D internal flows. An efficient ap- 
pr oach to solve fully consistent SOLA inversions was introduced 
by Jackiewicz et al.| (|2Qll ). In this paper we focus on inversions 
for three-dimensional vector flows on supergranular scales in the 
near- surface layers of the solar convection zone. 



1.1. Validating helioseismic inversions 

Two approaches have been used so far to validate time-distance 
inversions. The first approach consists of generating synthetic 
travel-time maps by convolving a (known) arbitrary frozen 
flow field with travel-time sensitivity kernels and then testing 
the inversion method using these synthetic travel times (e.g., 
Kosoviche v & Duvall| 1 1 997| ) . The second approach (Z hao et al. 



e.g., 



Gizon et al. \2009\ Moradi 



2007, 2010 ) is to use evolving realistic numerical simulations of 
three-dimensional radiative convection, where the helioseismic 
waves are naturally excited by the convection. 

The first approach is convenient, but may not represent a re- 
alistic situation, in particular because the prescribed subsurface 
structures are often too idealised. The second approach is prefer- 
able, but is limited by computing resources: only simulations 
of the very near- surface regions of the Sun are available today. 
Both approaches have shown that inversions are generally able 
to retrieve the horizontal components of velocity at supergran- 
ulation scales in the quiet Sun. However, some problems have 
been reported. For example, vertical velociti es have been mea- 
sured with the opposite sign near the surface ( Zhao et al. 2007| ) 
and the structure and sign of the flow s around sunspots vary with 
the inversion method used (see 
|etaLl|20TQ| ). 

The aim of this paper is to validate a particular implementa- 
tion of SOLA inversions for time-distance helioseismology us- 
ing a mixed approach: we take a snapshot from a large-box real- 
istic simulation of solar convection ( Ustyugov, 2008 ) to generate 
realistic travel-time maps by convolution with Born travel-time 
sensitivity kernels. A realistic noise component is added to the 
travel times. The travel-time maps are then inverted using a mul- 
tichannel SOLA inversion and compared with the known flows. 
This approach allows us to investigate various types of biases 
in the results of the inversion and to develop a robust proce- 
dure to minimise them. The most serious bias comes from the 
natural correlations among the components of flows induced by 
mass conservation that translate into a cross-talk between the 
components of the inverted flow (as we shall see, this is impor- 
tant when retrieving the vertical flow). Using realistic numerical 
simulations of convective flow velocities as input is very useful 
to set the acceptable level of random noise of the inverted flow 
velocities. The use of travel times with realistic noise properties 
is important in order to derive realistic estimates of the noise of 
the inverted flow velocities. 



2. Synthetic travel times 

To construct synt hetic travel times , we use a realistic hydrody- 
namic simulation (Ustyugov 2008 ) of the solar convection. The 
computational domain is a box 20 Mm in depth and 60 Mm in 
each horizontal direction. The simulation provides us with a rea- 
sonably realistic description of flows in the upper convection 



zone V 
nates are defined as 

X = (r,z) , 



, v^). Throughout this paper the spatial coordi- 



(1) 



where r is the horizontal position vector and z is the height. 

We choose a snapshot (t = 500 min, see Fig. 16 in the on- 
line supplement) of the above mentioned simulation for our in- 
version tests. Following |Du vail et al.| ( |199 7T), we consider diff'er- 
ent types of travel-time measurements between a surface point 
at position r and a concentric annulus or quadrants in order to 
measure travel times sensitive to flows in inward-outward (de- 
noted by "oi"), west-east ("we"), and north-south ("ns") direc- 
tions. Travel-time maps are denoted by T^(r), where the super- 
script a is an integer that uniquely refers to a particular combina- 
tion of choices in the data analysis: a type of geometry (oi, we, or 
ns), annulus radius (from 7.3 Mm to 29.2 Mm every 1.46 Mm), 
and a wave filter (here, ridge filters for one of /, pi, P2, P3, or p4 
modes). Thus, index a refers to one of M = 3x16x5 = 240 pos- 
sibilities. Additional information describing the measurement 
procedure in time-distance helioseismology is given by |Gizon 
|& Birch] ( [2005] ). 

Travel-time maps (travel times as functions of position vec- 
tor r) are generated by convolving the convection snapshot with 
sensitivity kernels according to 



= f K\/-r, 

Jo 



z)-v(r\z) dVdz + ^^(r). 



(2) 



where K'^ = (K^,K^, K^) is 3. vector travel-time sensitivity ker- 
nel (see, e.g., Fig.[T7]in the on-line supplement), v is the velocity 
vector of convecting flows, and the volume integral is taken over 
the Sun. The noise component is denoted by n^. The original 
horizontal size of the simulation box of 60 Mm was too small 
for our purpose. Since the simulated velocities are periodic in 
the horizontal directions, we copied the simulation box 10 x 10 
times. 

Kernels ar e computed usin g the s ingle- scattering Born ap- 
proximation (Birch & Gizon[ 2007j ) and depend on eigen- 
modes of a background 1-D standard solar model (fChrist ensen- 
[Dalsgaard et aLj[1996J). The kernels are invariant under horizon- 
tal translations; horizontal averages of these kernels for various 
oscillation modes are displayed in Fig.[T8]in the on-line supple- 
ment to display their sensitivity in depth. All kernels used in this 
study have sensitivities only in upper-most 10 Mm of the con- 
vection zone. 

Solar waves are excited by the action of turbulent convec- 
tion in the Sun and therefore travel times are inherently noisy. 
Following Gizon & Birch" ("2004), we use a realistic noise co- 
variance matrix (see Fig. 19 in the on-line supplement). 



Aab(r - /) = E[n\r)n\/)l 



(3) 



to generate realizations of the noise. 

Each travel-time map has 400 x 400 pixels with spatial 
sampling of 1.46 Mm, corresponding to the pixel size of 
a Michelson-Doppler-Imager full-disc image (MDI; Scherrer] 
et al.[|1995] ). In total, we generate 240 difl'erent travel-time maps 



(one for each index a, see, e.g., Fig.[T]). These travel-time maps 
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Fig. 1. Examples of synthetic travel-time maps used in this paper 
(top row) and observed travel times from SOHO/MDI (bottom 
row). Temporal length of observation is T = 6 hours. 



have spatial power spectra that are similar to observed travel- 
time power spectra based on MDI data, as demonstrated in 
Fig. [2] However, our synthetic travel times have weaker power at 
low wave-numbers due to the complete lack of signal from the 
simulated convective velocities below kRQ = 60 (only simulated 
noise is present). 

We check to ensure that point-to-annulus sensitivity ker- 
nels and noise covariance matrices obey expected symmetries. 
Furthermore, the sensitivity kernels and covariance matrices are 
forced to decay smoothly to zero towards the edge of the compu- 
tation box by multiplying with a smooth spatial function having 
zeros far from the region of interest. In addition, we make sure 
that the horizontal integral of each sensitivity kernel for is 
zero at all depths. This is a consequence of symmetries associ- 
ated with point-to-point kernels for v^, which are insensitive to 
the mean value of vertical velocity ( [Birch & Gizon|[2007] ). 



3. Improved SOLA inversions 

3.1. Real-space formulation 

The inversion method we use here is an improved version of 
the Subtractive Optimally Localised Averages (SOLA) method 
described by ,Jackiewicz et al.| ( |2007l|2QQ8|[2QTT] ). 

The SOLA algorithm describes how to optimally combine 
a given set of travel time measurements to infer the underlying 
properties of medium. In this paper, we wish to retrieve vjj^^(jco), 
i.e., an estimate of the a-component of the flow velocity in the 
neighbourhood of position xq = (ro,Zo) in the solar interior. In 
practice, we search for a linear combination of the travel-time 
measurements. 



(4) 



where w^(r/) are inversion weights to be determined. Combining 
Eqs. ^ and (|4]), we have 



ro,z;zo) • v(x) d rdz 



Jo 

+ Xi^^^^^""^oUo)^''(^/), 



(5) 




kBQ 
p4/ns/29.2 Mm 




Fig. 2. Comparison of azimuthally averaged spatial power spec- 
tra |f^(^)p of synthetic (solid) and observed (dashed) travel times 
displayed in Fig.[T]as a function of kRQ, where k is the horizontal 
wave-number and Rq is the radius of the Sun. 

where = ('K^,'Ky,'K^) is a vector averaging kernel with 
components 



(6) 



The expression (firdz means d^x. 

The expectation value of inverted velocity vJJ^^(JCo) may be 
expanded as 



E[vr] 



invOS) 



where 



,;invOS)/ 



=1 



(xo) ^ 7C(r - ro, z; zo)v^(x) d'r dz . 



(7) 



(8) 



and vjf^ = vjf^^^\ The form of equation ^ will be useful later in 
this paper to optimize the inversion. Notice that the second term 
on the right side of ^ represents leakage of other flow compo- 
nents into the inverted one. This cross -talk may be a significant 
source of bias and needs to be studied. 

We search for inversion weights so that the vector aver- 
aging kernel '7C^(r - ro,z;zo) resembles a user-supplied target 
T^{r - ro, z; Zo). Because we wish to invert for the a-component 
of the flow, we choose the target with components: 



To{r - ro, z; zo) = T(r - ro, z; zo)SaB 



(9) 
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Table 1. Notations. 



Symbol 


Meaning 


Formula 


Va(x) 


Qf-component of flow velocity at x (from the simulation) 


T{r - To, z; zo)Va{x) dV dz 




Qf-component of flow velocity at position Xq targeted by the inversion 




contribution of the /^-component of flow velocity to the inverted a- 


/^7g(r-ro,z;zo)v^(x) dVdz 




component of flow velocity 






inverted of-component of flow velocity (contains a noise component) 


2/,a <(n - ro; Zo)r%n) = 2^ v^^^ + noise 



where T(r-ro, z\ Zo) is a function that peaks around xq = (ro, zo). 
diminishes rapidly away from that point, and has unit integral. 
dap is a Kronecker 6. Throughout this paper, we choose a simple 
Gaussian 



(4 In 2)3/^ 
T{r, z\ Zo) = — ^^P 



41n2. 



41n2 



(10) 



where S}^ and are full widths at half maximum (FWHM) of the 
Gaussian in horizontal and vertical directions respectively. 

A successful inversion will return a value of vjj^^(jco) close to 
the target velocity 



vj^xo) = I T(r -r^,z\ Zo)Va(^) dz , 



(11) 



and cross-talk velocities vjf^^^(xo) close to zero for p a. 

A list of the notations used for the various flow velocities 
referred to in the inversion procedure is given in Table [T] 

The problem to be solved to obtain the belongs to the 
class of constrained regularised optimisations. The terms to be 
regularised are the following: 

- The misfit, i.e., how far the averaging kernel is from the de- 
sired target function, 

;rV"uo)= r [7C>,z;zo)-T(r,z;zo)]'dVdz, (12) 
J© 

where TC^ is an implicit function of the weights w^. 

- The cross -talk quantifying the leakage of the signal from 
other flow components into the inverted one, 

Xr(w";zo) = X r [^^"(^'^'^o)]' dVdz . (13) 

- The variance of vjf^, which corresponds to the root-mean- 
square of noise-related fluctuations in inferred velocities. We 
refer to this quantity as predicted error hereafter, 

^^^(w^; zo) = Yj Zo)Kb(ri - rj)wl{rj\ zo) . (14) 

i,j,a,b 

- The ad-hoc term quantifying the spread of the weights in 
space around ro. 



^(w";zo) = ^[w>,-uo)]'. 



(15) 



The regularisation based on S ensures that weights decrease 
towards the edge of the horizontal domain in order to enforce 
the spatial locality of the inversion and to prevent the weights 
from oscillating in the spatial domain. 



In practice, we search for the weights that minimise the cost 
function 



Jo 



[7C>,z;zo)-r(r,z;zo)r d^rdz + 
+ r [7C;(r,z;zo)]'dVdz + 6^[w>,;zo)f + 

P^a a,i 

+ ^«(^^'^o)^«^(^^" " 0)<(0'^o). (16) 

i,j,a,b 

Trade-ofl' parameters yu, v, and e balance these terms. The strat- 
egy to be empl oyed in order to set these values will be described 
in Section [33l 

The inversion is subject to constraints 



I 



dV dz 7C(r, z; Zo) = S^p for aU p , 



(17) 



in order to scale the amplitude of the inverted flow v^^^ appropri- 
ately. 

By taking derivative of ([16]) with constraint ([T7| added with 
respect to the weights, the problem can be cast into a linear in- 
verse problem, as explained by Pijp ers & Thompson| (|1992 ) or 
Jackiewicz et al.| ( [20111) . In real space, the matrix to be inverted 



has {N^M+Py ^ (100 elements, where N = 200 is the number 
of grid points in one horizontal direction, P = 3 is the number of 
physical unknowns (three velocity components), and M = 240 
as already defined earlier. 

3.2. Fourier-space (Multichannel) formulation 

The full problem written in real space is intractable to be solved 
using no wadays computers , becau se the matrix to be inverted is 
too large. Jackiewicz et al. (20lT) found a solution to this prob- 
lem by transforming to spatial Fourier space, where the inverse 
problem decouples as a consequence of the horizontal translation 
invariance of the sensitivity kernels. Thus, instead of a big linear 
inverse problem (inverting matrix having (10^)^ elements), we 
solve 40000 small linear inverse problems (inverting matrices 
having 240^ elements) in wave- vector space. This a pproach is 
called a multichannel inversion ( Jensen et al. |1998| ). Here we 
summarise results of |Jackiewicz et al. ( [2011 ) for the sake of 
completeness. 

Following Jackiewicz et al. ( 201 1), we use the following def- 
inition of the Fourier transform, such that any function f(r) and 
its 2-D spatial Fourier transform f(k) are related according to 



m = hlY,f(k)^xp(ik-r), 
k 

m = A^y/(r)exp(-/A:.r), 



(18) 
(19) 
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where hx = 1.46 Mm and hk = 0.022 rad/Mm are the grid spac- 
ings in the real and Fourier domains respectively. 

For each non- vanishing wave-vector k = (kx, ky), the vector 
of weights W(k) = [w^(A:) w^(A:) . . . w^(A:)]'^ is the solution to 
the matrix equation, 

hlN^ A(k)W(k) = T(k) for A: 9^ . (20) 
Each matrix A(k) has M x M elements A at given by 

Aabik) = (2nfj dzK*{k,z)Kl{k,z) + 

— CO 

+ 00 

+ iiKabik) + edab • (21) 

The vector T(k) = [ti(k) hik) . . . tM(k)]^ has M elements 4 
given by 



ta(k;zo) = (2: 



-too 



(k,z)T(k,z;zo). 



In addition, the matrix equation for the case A: = is 

hlN^A(0) C 


where C is an M x P matrix with elements 



[ W(0) 




hlrm 


[ L 




U^/(hlN^) 



c. 



Jq 



(22) 



(23) 



(24) 



L is a 1 X P vector of Lagrange multipliers, and the vector is 
a 1 X P unit vector with components = dap- 



Solutions to the A/^ - 1 equations ( 20 ) and the equation ( 23 ) 
give Fourier components of the weights, w^{k). By taking an in- 
verse Fourier transform we obtain w^(r). Equation ^ then gives 
an estimate of the a-component of the flow, vJJ^^. 

The diff'erence between the above equations and those 
in [Jackiewicz et al.| ( |2011| ) arises from additional constraint 
terms (|13|) and (|1 SjiT^ 

3.3. Picking trade-off parameters 

Trade-ofl' parameters yu, v, and e control the balance between var- 
ious terms in the cost- function ([16]). 

In practice, during the inversion for each flow component, 
we compute a grid of solutions by varying all three trade-ofl' pa- 
rameters. For each solution on this grid, we compute the misfit 
(12), predicted error of results (14), amount of cross-talk (13), 
and the spatial power of weightsTTSl. All four quantities can be 
computed in spatial Fourier domain. 

Standard optimisation me thods to achieve optimality, such 
as L-curve analysis (see, e.g., |Hanseii| |1998( [Jackiewicz et al. 



|2008 ), are not particularly useful to our problem. The elbow of 
the L-curve, which is considered an optimal point in parame- 
ter space, was located where the predicted noise level was very 
large. We therefore developed our own strategy of selecting val- 
ues for the trade-off' parameters. 

We start with parameter yu, which controls the trade-off' be- 
tween the misfit and random error. We select ji such that the 
noise level of the inverted quantity is less than some target value. 



Table 2. Predicted inversion errors and the expected magnitude 
of the velocities at the three diff'erent depths discussed through- 
out this paper. 



Depth [Mm] 


1 3.5 5.5 


FWHM Sh [Mm] 
FWHM s, [Mm] 


15 15 15 
1.1 2.2 3.5 


T - A days 


Noise (Tx [ms~^] 

Inverted signal {v^^^^)]^^ [ms~^] 

Targeted signal (v^^^ )]^^ [ms~^] 

Noise (Tz [ms~^] 

Inverted signal {Vl""^))!^ [ms-^] 

Targeted signal (v^^^ [ms~^] 


14 20 28 
35 20 13 
40 25 20 
3 13 133 
5 5 5 
5 7 7 


10^ X r = 6 hours 


Noise cr^ [ms~^] 

Inverted signal {v^^^^)]l^ [ms~^] 

Targeted signal (v^^^ [ms~^] 

Noise (Tz [ms~^] 

Inverted signal (v^"^^)^^ [ms-^] 

Targeted signal (v^^^ [ms~^] 


7 9 7 
41 25 19 
40 25 20 
2 2 1 
5 6 4 
5 7 7 



Notes. Presented results are for inversions with cross-talk minimised 
and for two discussed cases: (1) inversion using travel-time maps aver- 
aged over few days and (2) averaging over many realisations of similar 
flow structure, each averaged over a short time (here we assume 10"^ 
realisations, each averaged over 6 hours). 



chosen by the user, depending on the problem at hand. For ex- 
ample, one may choose a target noise level of 20 m/s to invert 
for horizontal flows in supergranules. A much lower noise level 
will be required to invert for the vertical component of the flow 
( [Duvall & Birch[ [20 1 Q| ) . The example noise levels cTq, are given 
in Table [2j where we give the predicted inversion error for v^^^ 
and v^^^ and the root-mean- square of the flow averaged with the 
target function (v^^^ and v^^^) and averaged with the resulting av- 
eraging kernel (yf^ and v^^^). It may be that this selection of jd 
leads to an averaging kernel which does not resemble the desired 
target function at all, in which case one would have to allow for 
a coarser spatial resolution. 

For the given yu, we then choose v so that the degree of cross- 
talk is less than 10"^ Mm~^ . We observe that this constraint 
places the following upper bound on the magnitude of the cross- 
talk terms: 



max 


Z 7C-(r,z;zo) 







max |7C^(r,z;zo)l 



< 0.05 . 



(25) 



The parameter e controls the degree of spatial confinement 
of the inversion weights and has an impact on the misfit and error 
as well. If 6 is too large, the weights will be highly localised in 
space around the central point, but a large misfit will result. 

For the particular cases that we study here, we find that local- 
isation of weights is accomplished when 5" < 2 x 10"^ km^s"^. 
The impact of this regularisation term is illustrated in Fig.|3] For 
the case shown, weights are confined to within a disc of radius 
-50 Mm. 
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Fig. 3. Example inversion weights w^(x,y) for 1 Mm depth and effect of the minimisation of the spread S . A strong regularisation e 
confines the inversion weights in the spatial domain (bottom row, e = 10^-^) compared to the case 6 = (top row). The cut aty = 
(horizontal black line) is displayed on the right-hand- side panel, where the values of w^(x,y = 0) were scaled by w^(x = 0,y = 0). 



3.4. The code 

The code implementing the above described procedure is writ- 
ten in Matlab. Matlab provides a compromise between compu- 
tational efficiency and the availability of higher-level software 
constructs that make the code lucid, modular, and easy to mod- 
ify. MATLAB-based code may be compiled into binary-executable 
form suitable for pipeline pseudo- automatic processing. The in- 
put is a text file, in which the user specifies all relevant parame- 
ters including the kernels to be used in the inversion and the set 
of trade-off' parameters to be investigated. This allows for the de- 
velopment of user interfaces such as the web-based or graphical 
interface, which can serve as the user front-end of the inversion 
code. 

The code is parallel and scales linearly with number of par- 
allel jobs. The execution is fast, e.g., a set of inversions for one 
ffow component involving 240 kernels, 200x200 spatial points, 
and a grid of 200 trade-off' parameter values takes around 6 
hours using 48 Opteron-2.3GHz CPUs. When many CPUs are 
involved, the extensive input-output load becomes a bottle-neck 
and aff'ects the total processing time. 

4. Inversion for horizontal flow 

We applied the above-described method to known synthetic 
travel times in order to validate the set-up and the performance 
of the inversion procedure. We focus ffrst on the horizontal ffow 
components, i.e., a = x or a = y. We compare two principal 
quantities: (1) the ffow map that is obtained by convolving the 
known velocity ffeld with the target function (vjf^), i.e., the best- 
case inversion scenario, and (2) the ffow map actually resulting 
from the inversion (vJJ^^). In this manner, we may investigate in 
detail diff'erent sources of bias in the results. 

The simulation convolved with the target function gives us 
from which we estimate the expected magnitude of the ffow 
that we want to invert for (see Table |2]). This places limits on 
the required noise level of the inversion so that the results have 
signal-to-noise ratios larger than 1 . Requirements on the targeted 
error level ffx the trade-off' parameter ji. 

The choice of the target function depends on the discretion 
and needs of the user. Here we focus on layers in the top few 
Mm, in particular on depths -zo = 1, 3.5, and 5.5 Mm. This set- 
up was selected because similar ffow inversions were also per- 



formed and discussed by |Jackiewicz et al.] ( |2008| ), which makes 
it possible to compare the results of both methods. As discussed 
in the preceding sections, the outcome of the inversion is a set 
of weights w^(r), which are used to combine the travel-time 
maps in order to obtain estimates of velocity V^^(r) and v^y^(r). 
Example weights for vf^ are displayed in Fig.js] 

All components of the averaging kernels 7C^(r,z;zo) are 



shown in Figs 20-22 (available in the electronic supplement). 
We also show comparisons between inversions when the cross- 
talk is minimised and not minimised. The minimisation of the 
cross-talk with v = 100 is very efficient. However, for depths 
of 3.5 and 5.5 Mm the minimisation of the cross-talk introduces 
some small artefacts in the inversion averaging kernel TCjJ, which 
are a small price to pay. 

Fig. |4] shows vertical cuts through the target functions and 
the averaging kernels when the cross-talk is minimised. Except 
for the target depth of 5.5 Mm, the depth dependences of the 
averaging kernels "K^ resemble that of the target functions. The 
near-surface inversion at 1 Mm is dominated by the /-mode. The 
inversion at 3.5 Mm has equal contributions from / and pi, with 
opposite signs, as shown in Fig. [23] in the on-line supplement, 
where we plot the contribution of individual modes to the aver- 
aging kernels "K^ as a function of depth. 



-x 10 




-10 -5 
z[Mm] 



Fig. 4. The cut through the x = y = point of the averaging 
kernel (solid) and the respective target function (dashed) for Vx 
inversion with minimised cross-talk, using travel times averaged 
over 4 days, at three discussed depths (1, 3.5, and 5.5 Mm). 
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The validation of the V;^ inversion is demonstrated in Fig. [5] 
Here we plot the desired v^^^ at three different depths and the in- 
verted vf^ without noise contributions. These two are very close 
for the depths 1 and 3.5 Mm. The differences between v^^^ and 
vf^ at depth 5.5 Mm are caused by an imperfect averaging ker- 
nel. When the random noise is added to the solution (bottom 
row of Fig. [5]), we see that the inversions for Vx at 1 Mm and 
3.5 Mm are still very good, while the inverted Vx at depth 5.5 Mm 
is dominated by noise. Minimising cross-talk does not improve 
the quality of the solution in this case. 

In order to quantify the inversion biases, we compare directly 
the expected and inverted values for a set of spatial locations. 
We select points separated by 7.5 Mm in each horizontal coor- 
dinate. This sampling interval is equal to half of the horizontal 
FWHM of the target function (thus the points are somewhat in- 
dependent). In addition, to avoid possible edge effects, we cut 
the outer part of the horizontal plane so that only the central 
200x200 pixel patch is kept. The scatter plots comparing vari- 
ous inversion components are displayed in Fig.|6]for the depth of 
1 Mm, the plots for the depth of 3.5 Mm are qualitatively similar. 
We estimate that for the depth of 1 Mm, the imperfect averaging 
kernel leads to an average underestimation of the horizontal flow 
components by some 20% (Fig.[6]left). Furthermore we note that 
in the case of the inversion for the horizontal flow components, 
the bias caused by the cross-talk is not important and the random 
noise level corresponds to the predicted value (Fig. [6] right). The 
results for the depth 3.5 Mm are similar, while the results for the 
depth 5.5 Mm are dominated by the random noise. 

Inversions for at the depth of 5.5 Mm are dominated by 
noise; this is evident from the plot of the azimuthally averaged 
power spectra of the random noise and of the signal as a func- 
tion of kRQ (Fig. |7] right). For comparison, similar plots for the 
depths 1 Mm and 3.5 Mm are given also in Fig. [7] where the 
power spectrum of the signal is well above the random noise 
near supergranular spatial scales. The decrease in power of the 
signal at low kRQ in Fig. [7] is because the convection simulation 
does not contain these scales. 

We showed that it is possible to retrieve Vx and v^; in the top 
3.5 Mm without noticeable bias and within the predicted noise 
level of ~25 m/s (for observing time T = 4 days) with 240 differ- 
ent travel-time measurements for ridges / to p4. In Table [3] we 
summarise our findings. We state the statistical quantities (cor- 
relation coefficient and the slope of the linear fit) comparing vJJ^^ 
and the corresponding vjf^. We also show the signal-to-noise ra- 
tio of the inverted vjf^. 

The cross-talk is unimportant in inversions for the horizontal 
flow components. The cross-talk could come only from v^, which 
is weak. We estimate that the influence of the cross-talk is less 
than 5%. The inversion at 5.5 Mm is already dominated by noise. 
If we drop the requirement on horizontal resolution, inversions 
at a depth of 5.5 Mm are also possible (e.g., with a FWHM of 
Sh = 25 Mm). 

It is possible to perform inversions with less temporal aver- 
aging, but to obtain a reasonable signal-to-noise ratio, we would 
have to relax the demand on the match between the averaging 
kernel and the target function. Typically, this leads to side-lobes 
in the averaging kernel (especially in the z-direction), which 
may make interpreting results more difficult. Another possibility 
would be to include a larger number of independent travel-time 
measurements (and therefore more sensitivity kernels). These is- 
sues are being worked on. 



5. Inversion for vertical flow 

5.1. Specificity of inversions for vertical flow 

Inversions for require a different methodology because sen- 
sitivity kernels have zero horizontal integrals at each depth 
(and therefore zero total integral) implying (see Eq.([6])) 



+ OC 

/ 



dV7(3(r,z';zo) = 0, foraUz'. 



(26) 



Consequently, with measurements discussed here, it is impos- 
sible to retrieve horizontal average (v^) of vertical flow and we 
may only invert for fluctuations - (v^). Equation ([26]) implies 



r dVdz'7CJ(r,z';zo) = OforyS: 
Jo 



(27) 



and thus equation ^Tl) cannot be written for a = z. This implies 
in turn that the matrix equation ( 23 ) for A: = may not be written 
either. 



5.2. Ignoring k = 

A solution to this issue is simply to replace Eq. ( [23] ) by 

<(0) = for aU a . (28) 

Since the averaging kernel must have zero integral, the av- 
eraging kernel will be offset by a small negative constant away 
from the central peak. Luckily, regularization by the term (15) 
will ensure that the averaging kernel will drop to zero at the 
edge of the inversion box. However, this solution is not quite 
satisfactory, because we do not have very much control over this 
extended negative surrounding sidelobe. 

5.3. Target function with zero mean 

It is more elegant to select a target function T{r - ro, z; Zo) with 
vanishing horizontal integral at each depth. If this target func- 

Table 3. Correlation analysis of various inversion components 
using synthetic travel times averaged over 4 days. 



Crosstalk ignored 



Depth [Mm] 


1 


3.5 


5.5 


corr(v^7,vj,"^) 


0.93 


0.75 


0.54 


slope(v7,vf ) 


1.03 


1.04 


1.17 


SNR{^^^') 


2.63 


1.10 


0.54 


corr(vj"\vj"^) 
slope(vj"\ vj"^) 
SNR{\^'^'') 


-0.63 
-1.98 
2.01 


-0.02 
-0.07 
0.42 


0.06 
1.67 
0.04 



Improved inversion, cross-talk minimised 



Depth [Mm] 


1 


3.5 


5.5 


corr(v^"\v^"^) 


0.92 


0.73 


0.51 


slope(v^7,v^7) 


1.02 


1.05 


1.29 


SNR{\^'^') 


2.40 


1.03 


0.46 


corr(v^"\v^"^) 


0.82 


0.31 


0.05 


slope(v^"\v^"^) 


0.92 


0.91 


1.41 


SNRi^l"") 


1.67 


0.35 


0.03 
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Fig. 5. Comparison of inverted with input 
data. Top-row panels show the input flow field 
convolved with the target function and the av- 
eraging kernel respectively. Bottom rows show 
inversion results in cases when cross-talk is ig- 
nored (left) and minimised (right). Random er- 
rors of the inversion are given in Table |2] 
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-101 
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Fig. 6. Inversion biases for at 1 Mm depth. 
Left: Noiseless v^"^ with cross-talk minimised 
versus the ideal vj\ The departure from slope 
unity (grey dashed line) is due to an imperfect 
match between TC and T. Right: Noisy v^/^ ver- 
sus noiseless v^"^. The results are plotted in two 
cases: when the cross-talk is ignored (grey x) 
and minimised (black +). The linear fit to the 
black crosses coincides with the dashed line of 
slope unity. The black dotted lines represent the 
predicted error of 14 ms~\ which is consistent 
with the observed scatter of the black crosses. 



tion is to peak around jcq = (ro, Zo), then this peak must be com- 
pensated by negative side-lobes in horizontal directions. These 
side-lobes have to be constructed in such a way that they do not 
lead to significant biases. We suggest to choose 



nr,z;zo) = H(r) 



■ exp 



41n2 



-(z-zoY 



where 



H(r) = 



41n2 



exp 



41n2. 



(29) 




c = exp - 



(30) 



(31) 



ensures that the horizontal integral of 7~ is zero. The horizontal 
part of the target function is constructed from a Gaussian peaked 
at ro = 0, from which a wide surrounding Gaussian annulus is 
subtracted. The free parameter n balances the width of a side- 
lobe compared to the width of the central peak. 
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Fig. 7. The azimuthally-averaged power spectra of various components of the Vx inversions using travel times averaged over 4 days. 
For reference, we plot the power spectrum of vf^ (thick grey solid line). We plot also the power spectrum of vf^ (solid line) and the 
power spectrum of the noise (i.e., the power spectrum of v^^^ - v^^^; dashed line) for the inversion where the cross-talk is minimised 
(black) and ignored (grey). 



The inversion is performed for each A: ?^ with this new 
target function, together with w^(A: = 0) = 0. 



5.4. Both target functions provide similar answers 

We compared the results obtained using the two target functions 
proposed above and found two conclusions. 

The two solutions described above provide results which 
are very close (see, e.g.. Fig. 24 in the electronic supplement). 
However, the solution with the target function having zero hor- 
izontal average is elegant and provides more control over the 
solution averaging kernel. 



By comparing the results obtained with various sizes of the 
negative annulus surrounding the central peak we found that the 
bias in the inverted flow caused by this negative sidelobe is neg- 
ligible if ^ > 3 in Eq. ( [30| . 

The results presented in the following sections are obtained 
using the first formalism. 






A X [Mm] 

Fig. 8. Horizontal spatial cross-correlation between V;^ and 
from the simulation in the top 1 Mm. These flow correlations to- 
gether with a non- vanishing cross-talk averaging kernel TC^ may 
lead to a systematic bias of the inferred v^. The circle of radius 
of 15 Mm shows the width of the desired target function. 




z[Mm] 



Fig. 9. The cut through the x = j = point of the averaging 
kernel (solid) and the respective target function (dashed) for 
inversion using travel times averaged over 4 days. Compare with 

Fig. a 

5.5. Validation of vertical flow inversion 

As is evident already, inversions for vertical flow are not as 
seamless as in the case of horizontal components. This is mostly 
because the vertical flow is much weaker on supergranular scales 
than the horizontal flow. As shown in Table [2j the expected RMS 
of the vertical flow in the top layers is of the order of 5 ms"^ 
Therefore, the predicted noise of the inverted has to be set to a 
much smaller value than for the horizontal velocities. As a result, 
the match between the desired target function and the resulting 
averaging kernel will have to be poorer. 

Furthermore, it is absolutely crucial to minimise the cross- 
talk to avoid the leakage of the large horizontal velocities into 
the small inferred vertical velocity. Even an apparently negligi- 
ble cross-talk averaging kernel could in the end cause a signif- 
icant bias in the results. Minimising the cross-talk is especially 
important because of the natural correlations between the verti- 
cal and horizontal flow components in the mass-conserving flow 
of the supergranules (see Fig. [8]). In the upper layers, horizontal 
outflows are associated with upflows. The structure of the cross- 
talk averaging kernel as shown in the top half of Fig. fTolin 



the case when cross-talk is not minimised then implies a negative 
bias in V, due to v^^^(^> and v^^^^\ 

Fig.[9]shows vertical cuts through the target functions and the 
averaging kernels when the cross-talk is minimised. The averag- 
ing kernels "K^ for the inversions at depth 3.5 and 5.5 Mm have 
sidelobes towards the surface. All components of the averaging 
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Fig. 10. All components of the averaging kernel for inversion at 1 Mm depth with a FWHM of = 1.1 Mm and Sh = 15 Mm. 
Bottom row: with cross-talk minimised, top row: cross-talk is ignored. The cross-talk is presented in the form of and TC^ 
averaging kernel components. Random error of the results is 3 ms"^ when assuming data averaged over 4 days. Over-plotted 
contours, which are also marked on the colour bar for reference, denote the following: half-maximum of the kernel (white), half- 
maximum of the target function (black), and by grey lines ±5% of the maximum value of the kernel (solid and dotted, respectively). 
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Fig. 11. Comparison of inverted with input 
data, compare with Fig.js] (v^) = 11 ms~^ 



kernel 7CJ(r,z;zo) for are displayed in Fig. 10 (and Figs. 25 
26| available in the electronic supplement), again comparing me 
cases when the cross-talk is and is not minimised. The action of 
the cross-talk minimisation term is very efficient. 

The validation of the inversion is demonstrated in Fig. 1 1 
Here we plot v/ at three different depths and the inverted v^^^ 
without noise contributions. These two are very similiar at the 
depth of 1 Mm. Small differences between v^^^ and v^^^ caused by 
an imperfect averaging kernel are visible at the depth of 3.5 Mm, 
become more significant at depth 5.5 Mm. When random noise 
is added to the solution (bottom row of Fig. [TT]), we see that 



the inversion for vertical flow is possible at 1 Mm depth only if 
cross-talk is minimised. Inversions at depths 3.5 and 5.5 Mm are 
buried in random noise, where minimising cross-talk does not 
help. The magnitude of cross-talk at a depth 1 Mm is much larger 
than that of the vertical flow; further, cross-talk is highly anti- 
correlated with the vertical flow, which makes the correlation 
coefficient between v^^^ and v^^^ close t o - 1 . Note that a si milarly 
high anti-correlation was measured by |Zhao et al. ( 2QQ7| ). 

All components v^^^^^ and the random noise component in 
the inversion for at 1 Mm depth are displayed in Fig. 12 We 
see in the top row that the leakage of the horizontal components 
^^inv(x) yinvO;)^ covcrs up Completely the weak signal of the 
vertical flow when the cross-talk is not minimised. 

Biases in the inversion at a depth of 1 Mm may be quan- 
tified by directly comparing expected and inverted values for a 



set of spatial locations in the horizontal plane (Fig. [13]), as in 
the case of the Vx inversion (Section |4]). The imperfect averaging 



kernel in this case does not cause any significant bias (Fig. 13 



left). The inversion is not sensitive to the horizontally aver- 
aged vertical flow (v^), which is therefore subtracted from v^^^ 
in the corresponding plot. The effect of minimising cross-talk is 
shown in Fig. 



13 



is anti-correlated with 



expected from me inversion due to the leakage of the horizontal 
mass -conserving flow components into the vertical one. 

Meaningful inversions for the vertical flow are possible in 
shallow near- subsurface layers of the convection zone only when 
the cross-talk between the vertical and the horizontal compo- 
nents is minimised. Vertical flow inversions on supergranular 
horizontal scales at depths greater than ~1 Mm performed us- 
ing / to p4 modes and travel time maps averaged over 4 days are 
dominated by random noise and the signal of the vertical flow 
may not therefore be inferred at all (see summary in Table [3]). 
This fact is also demonstrated in Fig. 14 where we plot power 
spectra of individual inversion components as a function of the 
spatial scale. At depths larger than 1 Mm, the power of the sig- 
nal is much less than the power of the random noise regardless 
of the spatial scale. 

6. Beating the noise: statistical averaging 

As demonstrated in previous sections, it is very difficult to have 
a meaningful inversion for solar flows on supergranular scales 
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Fig. 12. All components of the inversion at 1 Mm depth. Top row with the cross-talk ignored, bottom row with the cross-talk 
minimised. We demonstrate that if cross-talk is not addressed, horizontal components will leak into the inverted and cause a bias. 
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Fig. 13. Inversion biases for at 1 Mm depth, 
similar to Fig. [6] While the black line (with 
minimised cross-talk) almost coincides with 
the dashed line of slope unity, the grey one 
(cross-talk ignored) indicates a horrible bias. 
The predicted error is 3 m s~^ The inversion is 
not sensitive to the horizontally averaged hor- 
izontal flow (v^) = 11 ms~\ therefore it was 
subtracted from v^^^ in this comparison. 




Fig. 14. The azimuthally-averaged power spectra of various components of the inversions using travel times averaged over 4 days. 
Compare with Fig. [t] for Vx inversions. Notice that in the case of v^^^ containing the cross-talk at 1 Mm depth, the excess in power 
around kRQ ~ 150 is not caused by the random noise, but by the bias coming from the cross-talk (not plotted separately). In the 
other two cases the random noise is a main cause why these inversions are not possible. 



averaged over a short time, even in the shallow near- subsurface 
layers of the solar convection zone. The signal is overwhelmed 
by random noise. It is not feasible to average over longer times, 
because the expected lifetime of convection on supergranular 
scales is on the order of a day. By averaging over longer time, a 
significant portion of the scientifically useful information would 
be lost. 

It seems feasible to solve the issue by averaging over many 
realisations of similar flow structures each averaged over short 
time (therefore noisy). This concept was already used by other 



authors (e.g., [Duvall et al. 



2006|). As an example we note the 



possibility to measure the flows in many individual supergran- 
ules and to average over this sample in order to obtain th e typical 
flow structure in an "average supergranule", similar to Duv all &| 
Birch (2010). The predicted error in the results after averaging 
scales as 1 / ^^N, where N is the size of the ensemble of inde- 
pendent flow realisations. This scaling law allows us to relax the 
noise constraint on the inversion and therefore regularise more 
strongly about the misfit and cross-talk terms. See Table [2] for 
required random error of inversions results. 
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Fig. 15. Comparison at three depths of sta- 
tistical averages of velocity-vector over many 
flow realisations. The horizontal flow is dis- 
played by arrows while the vertical flow is 
color-coded. The reference arrow indicates 100 
ms~\ the random errors of the inversions are 
given in Table[2] Compare to Figs.[5]andpT] 



For instance, let us assume that we average over 10"^ inde- 
pendent inverted flows, each obtained by considering travel-time 
maps averaged over 6 hours. Individual flow maps are noisy 
and therefore contain little useful information about underlying 
flows. This selection allows us to relax the constraint on the es- 
timated error level by a factor of a 100 and therefore obtain a 
much better fit to the target function (see Fig. [27] in the elec- 
tronic supplement). With this set-up, the validation displayed 
in Fig. [T5| shows almost perfect correspondence between ex- 
pected and inverted v^^^ flows in the top 5.5 Mm of the 
convection zone. Inversion for at 5.5 Mm with a FWHM of 
Sh = 15 Mm is an exception, because the selected target function 
cannot be matched by the averaging kernel. Difl'erent inversions 
for Vz at this depth are possible with greater averaging, e.g., with 
Sh = 25 Mm. It is still crucial to minimise cross-talk in the case 
of Vz inversions to retrieve the correct answer. Power spectra of 
inversion components at all discussed depths are displayed in 
Fig.[28]in the electronic supplement. 



7. Conclusions 

We improved and validated an inversion algorithm based on a 
SOLA inversion approach. The formalism, algorithm and the 
code is universal and can, in principle, be used to invert for any 
quantity describing inhomogeneities in the solar plasma, pro- 
vided that the corresponding sensitivity kernels are available. 
The code is also ready to be used for application to real measure- 



ments. It will become part of the helioseismic pipeline running 
at German Data Center for SDO. We plan to use it to analyse all 
available SOHO/MDI and SDO/HMI data in order to routinely 
provide a tomographic image of the structure of the solar upper 
convection zone. 

The code is absolutely scalable allowing to include more 
independent measurements (and therefore more sensitivity ker- 
nels) in order to further refine the precision of the results. Thanks 
to the decoupling, the problem will still be solvable even using 
nowadays computers. 

We improved the inversions by introducing additional terms 
that allow to control and minimise some sources of bias in the 
results. Most importantly, we minimise the cross-talk between 
individual flow components, which is crucial especially for 
inversions. In principle, the formalism allows to minimise cross- 
talk between any selected quantities during the general inver- 
sion. The validation performed here silently assumes that the 
sensitivity kernels and noise covariance matrices are perfect. For 
the validation of the method and the code it is important if the 
sensitivity kernels and noise covariance matrices are solar-like. 
This allows us to study difl'erent sources of biases we may expect 
in the Sun. It is also very important that the sensitivity kernels 
used in the inversion contain all the details of the travel-time 
measurements, including the instrumental function afl'ecting the 
solar oscillation power spectrum. 

We found that by considering / to p4 frequency-averaged 
modes and supergranular spatial scales it is possible to perform 
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reliable and trustworthy flow inversions of the travel-time maps 
averaged over a few days in the top 3.5 Mm layer of the con- 
vection zone in the case of horizontal and Vy components, and 
in the top 1 Mm in the case of vertical component. Based on 
our experiment, we expect that using travel-time maps averaged 
over 4 days, it is possible to measure 3-D velocities as weak 
as 10 ms"^ at the surface and horizontal velocities having am- 
plitude 20 ms"^ in the top 5 Mm. We estimate, that using the 
travel-times averaged over 1 day, it should still be possible to 
measure all components of the supergranular flow at the surface 
and its horizontal components in the top 5 Mm. The cross-talk 
minimisation is crucial in order to measure the correct vertical 
velocity. Its presence in the results may exp lain the opposite si gn 
of the vertical flow inversion discovered by [Zhao et al. (2007 ). 

By considering many flow realisations and statistical averag- 
ing we might go deeper in the convection zone to learn about the 
horizontal flow components and perform vertical flow inversions 
for the top few Mm depths. 
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Fig. 17. The point- to- annulus sensitivity Icer- 
nels for three flow components computed for 
the /-mode, distance 28 Mm and outward- 
inward geometry. The white circle represents 
the location of the averaging annulus. 
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Fig. 18. Horizontal averages of sensitivity kernels may serve useful when estimating 
which depths are easier to target. The trend for each mode/ridge was obtained by taking 
averaging over all as within the given mode. 
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Fig. 19. An example noise-covariance matrix for /-mode travel times averaged over 6 hours. 
In this plot, a stands for the combination of the /-mode, oi geometry, and annulus radius of 



[Mm] 7 3 Mm, b stands for the combination of /-mode, we geometry, and annulus radius of 8.8 Mm. 
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Fig. 20. All components of the averaging kernel for V;^ inversion at 1 Mm depth with a FWHM of = 1.1 Mm and su = \5 Mm. 
Bottom row: with cross-talk minimised, top row: cross-talk is ignored. Random error of the results is 14 ms"^ when assuming 
data averaged over 4 days. Over-plotted contours, which are also marked on the colour bar for reference, denote the following: 
half-maximum of the kernel (white), half-maximum of the target function (red), and ±5% of the maximum value of the kernel (blue 
and green, respectively). 
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Fig. 21. All components of the averaging kernel for v^^ inversion at 3.5 Mm depth with a FWHM of = 2.2 Mm and Sh = 15 Mm. 
Random error of the results is 20 m s"^ when assuming data averaged over 4 days. For details see Fig. 20 
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Fig. 22. All components of the averaging kernel for Vx inversion at 5.5 Mm depth with a FWHM of = 3.5 Mm and Sh = 15 Mm. 
Random error of the results is 28 m s"^ when assuming data averaged over 4 days. For details see Fig. 20 
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Fig. 23. The contributions of particular modes 
to the horizontally averaged averaging kernel 
for inversions using travel times averaged 
over 4 days for depths 1 and 3.5 Mm. We 
do not display the inversion for the depth of 
5.5 Mm, because it is heavily dominated by 
noise. 
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Fig. 24. To solve the peculiarity of the inversion, we introduced two formalisms in 
Section [5J] Here we plot performance of those. In black, the magnified section a j = and 
z - Zo of diff'erent target functions are displayed, one with removed mean (1) and one con- 
structed with negative side-lobes (2). The resulting averaging kernels are also plotted. It is 
evident that the resulting averaging kernels are qualitatively very similar even when diff'erent 
formalisms were used to compute them. 
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Fig. 25. All components of the averaging kernel for inversion at 3.5 Mm depth with a FWHM of = 2.2 Mm and Sh = 15 Mm. 
Random error of the results is 13 m s"^ when assuming data averaged over 4 days. For details see Fig. 20 
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Fig. 26. All components of the averaging kernel for inversion at 5.5 Mm depth with a FWHM of 
Random error of the results is 133 m s"^ when assuming data averaged over 4 days. For details see Fig. [20 



3.5 Mm and Sh = 15 Mm. 
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Fig. 27. The cut through the x = j = point of 
the averaging kernel (solid) and the respective 
target function (dashed) for the (left) and 
(right) inversions using averaging over many 
flow realisations plotted along with the cor- 
responding target functions at three discussed 
depths (1 Mm in blue, 3.5 Mm in green, and 
5.5 Mm in red). Compare to Figs. |4] and [9] 
where the resemblance of the target functions 
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Fig. 28. The azimuthally-averaged power spectra of the inversion components at depths of 1, 3.5 and 5.5 Mm for averaging over 
many flow representations. For reference, we plot the power spectrum of v^^^ using the black solid line. Then we plot the power 
spectrum of v^^^ (solid line) and power spectrum of the noise (i.e., the power spectrum of v^^^ - v^^^; dashed line) for the inversion 
where the cross-talk is minimised (blue) and ignored (red). Compare to Fig. [14] 



